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1. Introduction 

Wilson fermions [1] are a popular way of discretizing fermions in lattice QCD. Their main 
drawback is that they break chiral symmetry. As a consequence, they are afflicted with dis¬ 
cretization errors linear in the lattice spacing a. These 0(a) errors can be removed by adding 
the Sheikholeslami-Wohlert term (clover term) to the Wilson action [2] and further improvement 
terms to the matrix elements of interest. This method of 0(a) improvement has become known 
as ‘Symanzik improvement program’ [3, 4]. Moreover, currents which are conserved in a chiral 
theory and hence need not be renormalized require a finite renormalization with Wilson fermions. 
One of these currents is the isovector axial current which we are concerned with in this work. 

The isovector axial current is a quark bilinear and in its bare form it can be written as 

=FW7Ai75y (1-1) 

where T^‘ is a matrix acting in flavor space. Symanzik’s effective theory predicts that the axial 
current will mix with the derivative of the pseudoscalar density P^'(x) at 0(a), when the lattice 
discretization breaks chiral symmetry. This can be compensated for by adding a corresponding im¬ 
provement term to the bare current. Its coefficient is denoted by ca and is at the heart of improving 
the axial current. 

Furthermore, the axial current is renormalized by multiplying it by the renormalization fac¬ 
tor Za and another mass-dependent term. Thus, the fully renormalized and improved axial current 
on the lattice is 

(Ar)^(x) =ZA(l+bAamq) [Af^(x)+acAdfiP''(x)] with /^(x) = l/7(x) 75 y i/r(x). (1.2) 

The axial current plays a fundamental role in many applications, notably the computation 
of quark masses and meson decay constants in the pseudoscalar sector. These are not only of 
phenomenological interest, but they provide a precise way of setting a physical scale in lattice 
simulations, too. One observable that can be used for this purpose is the kaon decay constant /k [5]. 
In these contexts, it is crucial to employ the improved and renormalized current, since, typically, 
improvement and renormalization each contribute about 10-20% to the final result [ 6 , 7, 8 ]. Both 
can be computed in perturbation theory. However, in previous works it was found that the non- 
perturbative results deviate strongly from the 1-loop estimates. The deviations can be several times 
the 1-loop contribution itself. Therefore, it is desirable to determine ca and Za in a non-perturbative 
way. 

The methods we use to determine ca and Za have been introduced in previous papers, which 
applied them to the quenched [9, 10] and two-flavor QCD cases [ 6 , 8 ]. Here, we look at an action 
with 0(a) improved Af = 3 mass-degenerate dynamical Wilson fermions [11] and the tree-level 
Symanzik-improved gauge action (TLI gauge action) [12]. The main part of this text is about the 
renormalization of the axial current. The renormalization condition is summarized in section 2, the 
ensembles of gauge field configurations that we used are described in section 3, and our preliminary 
results, in particular for the interpolating function for Za (valid for lattice spacings < 0.09fm), can 
be found in section 4. The determination of the improvement coefficient ca was recently finished 
and will be published soon [13]. We will only reproduce the main result in section 4. A preliminary 
report can also be found in [14]. 
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2. Renormalization Condition 

Renormalization conditions for the axial current are based on the idea that they can restore 
chiral Ward identities, which are broken by the Wilson term, up to 0{a^). This is done by choosing 
one particular Ward identity and adjusting Za so that it holds exactly. The condition that we choose 
has been introduced for simulations with two dynamical fermions [8], but a similar condition has 
already been used in the quenched case [10]. We give a shortened description of its derivation 
below. More details can be found in the original papers. 

The Ward identity our renormalization condition is based on is similar to the PCAC relation, 
i.e., it is derived from a chiral rotation of the quark fields, but, in addition, the axial current A^(y) 
is inserted as an internal operator. The resultant identity is 

da^ (x) (a"^ {x)aI (y) ^ext^ - 2m d'^x (p‘^ (x)A^ (y) {V^ (y) ^ext), (2.1) 

where /? is an arbitrary region containing y, ^ext is an operator built from fields oufside /?, and Vf 
is fhe isovecfor vecfor currenf. As region /? we choose fhe spacefime volume befween fwo space¬ 
like hyperplanes. Furfhermore, eq. (2.1) can be modified by using fhe PCAC relafion fo shiff fhe 
infegrafion domains. Seffing v = 0 and confracfing fhe flavor indices a and b wifh one arrives 
af 


|d3xd3y£^^^(Ag(x)A(;(y)^ext) 

-Im jd^xd^y j dxo (^/^(x)Ao(y) ^ext^ 

= /|dV (V'o"(j)^ext) (2.2) 

wifh xq > yo defining fhe hyperplanes. 

We evaluate fhis idenfify on a laffice wifh Schrddinger functional boundary conditions (peri¬ 
odic in space, Dirichlef in time) [15, 16] wifh vanishing background field. The source operator ^ext 
is builf from fhe quark fields ^ and ^' af fhe boundaries xq = 0 and xq = T\ 

= ^23) 

6L° 

with 

^‘■ = a'^LC(u)75y«(u-v)C(v) and 75 y «(u- v) C'(v), (2.4) 

where fhe wavefuncfion ft) is undersfood fo be an approximafion of fhe pseudoscalar ground sfafe. 
Ifs consfrucfion is defailed in [13, 14]. We summarize fhe resulfs in section 4. The free index c 
fhaf appears in eq. (2.3) is confracfed wifh fhe free index from eq. (2.2). In fhis case, fhe ferm 
on fhe righf-hand side involving fhe vecfor currenf can be simplified fo fhe boundary-fo-boundary 
correlafor 

= (2.5) 
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up to 0{a^), as was shown in [10, 8] by using isospin symmetry. 

In the remaining two terms on the left-hand side, the continuum currents are replaced by their 
improved and renormalized counterparts from the lattice. The result can be put in this form, 

Z\ (l +Z7Aamq)^ [Fl^{xo,yo)-2m-Fp^{xo,yo)] =Fi, (2.6) 


with the improved correlation functions 

^AA (^ 0 , 3 ^ 0 ) = FAA{xo,yo) +aCA 4o^PA (^0,3^0) + 4 o^ap(^o,3'o) 
+ a^c\ djca dy^Fpp (ao , yo), 


(2.7) 


and 


^0 


^PA(-*o,yo) =a w(Ao)[fpA(Ao,yo)+«CA5vo7Tp(^o,yo)], 


( 2 . 8 ) 


f)=yo 


where d denotes the central difference operator and FxY{xo,yo) with X,Y £ {Ao,P} stands for 


FxY{xo,yo) = 


6L6 




(2.9) 


and 


w(ao) = 


1/2 if Aq = yo orxo = Ao 

1 if yo < ^ < .^0 


( 2 . 10 ) 


implements the trapezoidal rule. Note that in eq. (2.6) the renormalization factors, which would 
arise from the boundary fields ^ and cancel on both sides and that the product mP^' can be renor¬ 
malized with the same factor Za as the axial current due to the PCAC relation. The mass-dependent 
term proportional to bA will be dropped from here on, since we will impose the renormalization 
condition at vanishing mass. Of course, we can not tune the parameter to get exactly zero mass, 
but this only amounts to an 0{am) effect. Thus, our final renormalization condition is 


Za = lim 

m^O 


■_ F_x _■ 

.T’aaC-^oAo) -2m-Fj5^(xo,yo)_ 


1 

2 


( 2 . 11 ) 


In order fo maximize fhe disfance befween fhe inserfion poinfs, we choose xq = |r and yo = \T. 

Excepf for F\, eq. (2.11) is builf from correlators of fhe form given in eq. (2.9). When one 
evaluafes fhese by performing fhe Wick confracfions, one finds fhaf only six confracfions confribufe, 
which are illusfrafed in figure 1. Two of fhem are disconnecfed. As argued in [8], fhey only give 
rise fo O(a^) confribufions and cancel in fhe continuum limif. By omiffing fhem and faking only 
fhe connected confracfions, one obfains an alfernative renormalizafion condifion. We will denote 
fhe corresponding renormalizafion facfor by Z™". 
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Figure 1: The six non-vanishing Wick contractions of correlation functions Fxy with source and two 
bulk insertions X and Y, see eq. (2.9), taken from [8]. 


xr/a^ 


K 

#REP 

#MDU 

ID 

12^ X 17 

3.3 

0.13652 

10 

10240 

Alkl 



0.13660 

10 

13048 

Alk2 

163 X 23 

3.512 

0.13700 

2 

20480 

Blkl 



0.13703 

1 

8192 

Blk2 



0.13710 

3 

24560 

Blk3 

203 X 29 

3.676 

0.13700 

4 

15232 

Clk2 



0.13719 

4 

15472 

Clk3 

243 X 35 

3.810 

0.13712 

7 

15448 

Dlkl 


Table 1: Overview of simulation parameters, number of replica and total number of molecular dynamics 
units of gauge configuration ensembles labeled by TD’. 


3. Simulations 

To determine the renormalization factor Za of the axial current, we re-use the ensembles that 
were generated to obtain the improvement coefficient ca-^ These had already been designed to fit 
this purpose, in particular the ratio of the spatial and temporal extents was approximately T/Lta 
3/2 with L 1.2fm. From [18], we expect this to be a good trade-off between a large infrared 
cutoff and small O(a^) effects. An overview of the simulation parameters is given in table 1. 

The coupling j3 was chosen such that the physical lattice size L stays roughly constant (line of 
constant physics). In this way, O(a^) ambiguities in Za are guaranteed to smoothly vanish in the 
continuum limit. This was done using the perturbative relation between the lattice spacing a and 
the bare coupling. However, only the first two, universal coefficients bo and b{ of the beta-function 
could be taken into account, because higher-order terms are not known for the TLI gauge action. To 
test for deviations from the line of constant physics, the gradient (or Wilson) flow coupling gcF was 
computed [19]. It is a renormalized coupling that depends on L as a scale, i.e., it will be constant if 
L is constant. 

The parameter K was tuned towards a vanishing PCAC mass. In [11], an upper bound of 

*The ensembles were generated using the openQCD code [17], see also http : //luscher . web . cern . ch/ 
luscher/openQCD/. 
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lampcAcI < 0.015 was employed. We expect that Za is more sensitive to the mass. Since at most 
values of j3 we have several ensembles with different fc values, we can check its influence explicitly 
by comparing the results. 

In order to monitor the autocorrelation of the generated gauge configurations, several observ¬ 
ables were computed alongside which are defined in terms of the gauge field smoothed by means of 
the gradient flow [20,21]. Mostly, they showed autocorrelation times that did not exceed 250 MDU, 
only the topological charge became frozen at the largest j3 = 3.810. For ca and Za, this should 
amount to a cutoff effect, but we have estimated it explicitly for ca and have indeed found no 
significant deviations [13]. Further details, in particular about the algorithmic details can be also 
found in that reference. 

4. Results 

We have measured the correlators that are necessary to compute the renormalization factor 
and the PCAC mass on every second trajectory (in Alk2, A2kl) or every fourth trajectory (in the 
remaining ensembles). Via eq. (2.11), we compute Za as well as the alternative Z“", where only 
the connected Wick contractions are included. 

As already anticipated in eq. (2.4), our Schrddinger functional correlators involve boundary 
operators with a paiticular choice of wavefunction, which is constructed in such a way that it 
suppresses the contribution of the first excited state in the pseudoscalar channel. This optimal 
wavefunction was determined in the context of our non-perturbative calculation of the improvement 
coefficient ca in [13] (see also [14] for a preliminary report), which employs the same gauge field 
ensembles at constant physics and the same kinematical setup as used here. It relies on demanding 
the quark mass extracted from the PCAC Ward identity to stay unchanged when the external states 
are varied, where in practice these external states are modeled as superpositions of spatial trial 
(hydrogen-like) wavefunctions designed to approximately maximize the overlap with the ground 
and first excited state, respectively. For the purpose of Za, however, only the wavefunction for the 
approximate pseudoscalar ground state is required, which we therefore choose as the one already 
obtained in [13], i.e.. 


(O^m = '£rij°^(Oi, T](°) = (0.5317,0.5977,0.6000), (4.1) 

1=1 

in terms of the basis of suitable trial wavefunctions CO, mentioned above. With this wave function 
at hand, the eq. (2.11) for Zp^ is evaluated upon prior projection of all entering correlation functions 
to this approximate ground state. 

Moreover, as it is evident from eqs. (2.7) and (2.8), a genuine non-perturbative determination 
of Za also requires the knowledge of non-perturbative values for ca for our simulation parameters. 
We thus rely on the result of our aforementioned recent non-perturbative computation of ca in 
three-flavor lattice QCD with tree-level improved gauge action [13], which we reproduce here for 
convenience: 


CA(go) = -0.006033ggx [1+exp(9.2056-13.9847( 4 . 2 ) 


6 



Improvement and renormalization of the axial current in Ni = 3 LQCD 


Christian Wittemeier 



Figure 2: Final results for ca together with interpolation [13]. The dotted line indicates the perturbative 
1-loop asymptotics [22]. 


this formula is valid for bare couplings below gg « 1.8 and with statistical errors between ss 4% 
near the largest and ss 8% near the smallest bare couplings simulated. For further details about this 
determination of ca we refer to refs. [13, 14]. 

Our preliminary results for Za obtained for the ensembles considered so far are collected in 
table 2. In our tentative error analysis, we estimate the errors on Za for each replicum via a full 
autocorrelation analysis as described in [23] and compute a weighted average over all replica within 
an ensemble. 

Table 2 also includes the unrenormalized PCAC mass ampcAC, which is computed using the 
perturbative value of ca and the wavefunction ft)^(o), the gradient-flow coupling and the results 
for Z“" obtained via the alternative definition of the renormalization factor, which includes only 
connected contractions. 

As can be seen from the table, gQp is approximately constant, only the deviation on the ensem¬ 
ble Dlkl at the largest j3 is more pronounced. We do not expect that this deviation from constant 
physics has a significant effect on our result, but we plan to check it explicitly using the ensemble 
B2kl from [13], whose parameters are identical to the ones of Blkl except for j3. Some of the 
ensembles in the B and C groups show a significant though not yet severe mass dependence. We 
consider taking a closer look at this issue and adding new ensembles with different K values, too. 

For Lja = 12 (j3 = 3.3), Z“" differs significantly from the standard value Za, which seems 
to signal significant O(a^) uncertainties in Za at this lattice spacing {a k, 0.09fm). Similarly large 
cutoff effects at this jS were also observed for this action in [24]. However, at smaller lattice 
spacings the results for both definitions of Za are in good agreement within their errors. We plan 
to examine the impact of cutoff effects on Za more closely by adding an ensemble at L/a = 14 
(j8 = 3.414). 
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ID 

^JWIPCAC 

Sgf 

ycon 

Za 


Alkl 

-0.0010(7) 

18.12(21) 

0.8162(78) 

0.6553(90) 

* 

Alk2 

-0.0086(6) 

16.95(13) 

0.8290(92) 

0.6489(72) 


Blkl 

+0.0063(2) 

16.49(13) 

0.7757(23) 

0.7666(47) 


Blk2 

+0.0056(3) 

16.85(20) 

0.7758(45) 

0.7677(71) 


Blk3 

+0.0022(2) 

16.11(14) 

0.7804(30) 

0.7516(36) 

* 

Clk2 

+0.0066(2) 

15.53(14) 

0.7889(16) 

0.7888(48) 


Clk3 

-0.0005(1) 

14.64(13) 

0.7822(27) 

0.7785(33) 

* 

Dlkl 

-0.00269(8) 

13.90(11) 

0.7969(16) 

0.7904(16) 

* 


Table 2: Summary of results: the unrenormalized PCAC quark mass, the gradient-flow coupling, the results 
for Za using the alternative definition (Z™”) and the standard definition including disconnected contractions. 
Ensembles marked by are used in the fit procedure. 



Figure 3: Plot of our preliminary estimate of Za versus the bare coupling gQ. The data points and a Pade fit 
are shown. 


In figure 3, the Za values of the four ensembles with the smallest absolute PCAC mass at 
each value of jS are plotted against the bare coupling. They have been used to determine a Fade 
approximation of Z^igg) based on the ansatz 


ZA(g^o) = 


1+Qi -go+ai-go 
1 +^ 1 -go 


(4.3) 


which is constrained to yield the correct continuum limit Za( 0) = 1. The coefficients we found are 


= -0.6492, a2 = 0.0619, bi = -0.5298. (4.4) 
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With these parameters, the fit function lies well within the errors of all data points. 

Let us emphasize again that the results on Za should be regarded as preliminary, since a final 
error analysis as well as the inclusion of the reweighting factors (to compensate for the approx¬ 
imation errors of the RHMC algorithm employed for the third quark in our simulations) are still 
missing. Regarding the latter, our experience from ca leads us to expect only a minor influence 
from them. Moreover, to account for the topology freezing observed at the finest lattice spacing, 
we will supplement our definition of Za with the condition to restrict the analysis to the sector of 
zero topological charge, as we did for c\ [13]. 

5. Conclusions 

We have determined a preliminary expression for the renormalization factor of the 

isovector axial current for the tree-level improved gauge action and three dynamical flavors of 
Wilson fermions. It is summarized in the interpolation formula of eq. (4.3). Together with the 
improvement coefficient ca, which has already been determined non-perturbatively [13], this will 
make it possible to obtain precise results for matrix elements such as pseudoscalar decay constants. 
However, before our result is applied, we will scrutinize our analysis by projection to the zero- 
topology sector and also investigate deviations from the constant physics condition. In addition, a 
simulation at another point (L/a = 14, j8 = 3.414) along our line of constant physics is under way 
to shed light on the cutoff effects reflected at the coarsest lattice spacing by the differences between 
Za and Z™". 
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